三次埃尔米特插值


Ⅰ. 前言 –> Ⅱ. 作用 –> Ⅲ. 两种典型插值原理 –> Ⅳ. Python 代码


Ⅰ. 前言

涉及到部分拉格朗日插值基函数和牛顿插值均差系数，建议先阅读：

拉格朗日插值 ______ 牛顿插值


Ⅱ. 作用

埃尔米特（Hermite）插值是一类概念，可以认为不是某种具体的插值方法. 前面拉格朗日插值和牛顿插值只要求插值节点函数值与实际因变量相等即可，但在工程领域，像无人驾驶汽车的拐弯，这绝不是只要求函数值相等就可以了，如果曲率或者导数、导导数（二阶导数）甚至更高阶导数也能缓一点，那么转弯相对就平稳.

埃尔米特插值是要求插值函数与本源函数，在函数值角度相等，而且也在特点节点导数（或高阶导数）值也相等的一类概念. 为了避免龙格现象，工程上一般常用三次埃尔米特插值法已经能达到很好的效果.


Ⅲ. 两种典型插值原理

  1. 三点三次埃尔米特插值

自变量： (x_0 < x_1 < x_2)

因变量： (y_0,\ y_1,\ y_2)

满足条件：

  1. H(x_i)=y_i,\ \forall i\in\{0,1,2\}

  2. H'(x_1)=y_1'

定义插值多项式【除三次项系数被修正，其余符合牛顿插值 f[x_0,x_1,\ldots,x_n] 多项式形式】：

\displaystyle H(x)=f[x_0]\times 1+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+C(x-x_0)(x-x_1)(x-x_2)

由于在中间节点处导数相等，所以对 H(x) 求导后代入 (b.) 条件可得：

C=\displaystyle\dfrac{y_1'-f[x_0,x_1]\times 1-f[x_0,x_1,x_2](x_1-x_0)}{(x_1-x_0)(x_1-x_2)}

余项【其中 \xi=\xi(x)\in(x_0,x_2)】：

R(x)=y(x)-H(x)=\dfrac{f^{(4)}(\xi)}{4!}(x-x_0)(x-x_1)(x-x_1)(x-x_2)

  1. 两点三次埃尔米特插值

自变量： (x_0 < x_1)

因变量： (y_0,\ y_1)

满足条件：

  1. H(x_0)=y_0,\ H(x_1)=y_1

  2. H'(x_0)=y_0',\ H'(x_1)=y_1'

定义插值多项式：

H(x)=y_0\alpha_0(x)+y_1\alpha_1(x)+y_0'\beta_0(x)+y_1'\beta_1(x)

其中 \alpha(x) 和 \beta(x) 都是三次埃尔米特多项式基底, 且满足一定约束条件.

由于条件 (a.)(b.) 可以确定唯一的三次埃尔米特插值多项式（笔者只能从结果反推这个结论应该是正确的，严格证明找相应文献）. 所以有两种做法，一种是直接把条件代入三次多项式通式方程求解待定系数，还有一种是借鉴拉格朗日插值基底的性质.

\begin{align} l_{t}(x)=\dfrac{\prod\limits_{i=1}^{n\setminus t}(x-x_{i})}{\prod\limits_{i=1}^{n\setminus t}(x_t-x_{i})} \Longrightarrow l_t(x_s)= \begin{cases} 1,\ \ t=s \\ 0,\ \ t\neq s \end{cases} \notag \end{align}

基底满足条件：

\alpha_0(x_0)=1,\ \alpha_0(x_1)=0,\ \alpha_0'(x_0)=\alpha_0'(x_1)=0

\alpha_1(x_0)=0,\ \alpha_1(x_1)=1,\ \alpha_1'(x_0)=\alpha_1'(x_1)=0

\beta_0(x_0)=\beta_0(x_1)=0,\ \beta_0'(x_0)=1,\ \beta_0'(x_1)=0

\beta_1(x_0)=\beta_1(x_1)=0,\ \beta_1'(x_0)=0,\ \beta_1'(x_1)=1

接下来借鉴性（很有难度）地构造出基底范式：

\aleph(x)=(ax+b)\times\left(\dfrac{x-x_i}{x_t-x_i}\right)^2

通过范式模板，用待定系数把基底满足条件代入后可解出基底：

\alpha_0(x)=\left(1+2\times\dfrac{x-x_0}{x_1-x_0}\right)\left(\dfrac{x-x_1}{x_0-x_1}\right)^2

\alpha_1(x)=\left(1+2\times\dfrac{x-x_1}{x_0-x_1}\right)\left(\dfrac{x-x_0}{x_1-x_0}\right)^2

\beta_0(x)=(x-x_0)\left(\dfrac{x-x_1}{x_0-x_1}\right)^2 \ \ \ \ \beta_1(x)=(x-x_1)\left(\dfrac{x-x_0}{x_1-x_0}\right)^2

余项：

R(x)=\dfrac{f^{(4)}(\xi)}{4!}(x-x_0)^2(x-x_1)^2


Ⅳ. 代码

三点三次埃尔米特插值.py
 1'''
 2# System --> Linux & Python3.8.0
 3# File ----> 三点三次埃尔米特插值.py
 4# Author --> Illusionna
 5# Create --> 2024/2/20 21:18:17
 6'''
 7# -*- Encoding: UTF-8 -*-
 8
 9
10class HERMITE:
11    def __init__(self, X:list, Y:list, middlePointDerivative:float) -> None:
12        self.__X = X
13        self.__Y = Y
14        self.__middlePointDerivative = middlePointDerivative
15        self.__pos = False
16
17    def Coefficients(self) -> list:
18        """公有函数: 计算埃尔米特多项式系数."""
19        L = [self.__Y[0]]
20        tmpA = (self.__Y[1] - self.__Y[0]) / (self.__X[1] - self.__X[0])
21        L.append(tmpA)
22        tmpB = ((self.__Y[2] - self.__Y[1]) / (self.__X[2] - self.__X[1]) - tmpA) / (self.__X[2] - self.__X[0])
23        L.append(tmpB)
24        C = (self.__middlePointDerivative - tmpA - tmpB*(self.__X[1] - self.__X[0])) / ((self.__X[1] - self.__X[0]) * (self.__X[1] - self.__X[2]))
25        L.append(C)
26        self.__coefficients = L
27        self.__pos = True
28        return L
29
30    def Interpolate(self, x:float) -> float:
31        """公有函数: 插值."""
32        if self.__pos == True:
33            result = 1*self.__coefficients[0]
34            tmp = 1
35            for index in range(0, len(self.__X), 1):
36                value = x - self.__X[index]
37                tmp = tmp * value
38                result = result + tmp * self.__coefficients[index+1]
39            return result
40        else:
41            self.Coefficients()
42            return self.Interpolate(x)
43
44
45if __name__ == '__main__':
46    print('\033[H\033[J')
47
48    X = [1/4, 1, 9/4]
49    Y = [1/8, 1, 27/8]
50    middlePointDerivative = 1.5     # 中间点导数, 即 x1 处的导数.
51
52    obj = HERMITE(X, Y, middlePointDerivative)
53    print(f'埃尔米特多项式系数:\n{obj.Coefficients()}')
54    print(f'x = 1.6 处的插值结果: {obj.Interpolate(1.6)}')
两点三次埃尔米特插值.py
 1'''
 2# System --> Linux & Python3.8.0
 3# File ----> 两点三次埃尔米特插值.py
 4# Author --> Illusionna
 5# Create --> 2024/2/20 21:46:23
 6'''
 7# -*- Encoding: UTF-8 -*-
 8
 9
10class HERMITE:
11    def __init__(
12            self,
13            X:list,
14            Y:list,
15            derivative:list
16    ) -> None:
17        self.__X = X
18        self.__Y = Y
19        self.__derivative = derivative
20
21    def Interpolate(self, x:float) -> 'function':
22        """
23        公有函数: 返回插值函数地址.
24        """
25        (x0, x1) = (self.__X[0], self.__X[1])
26        (y0, y1) = (self.__Y[0], self.__Y[1])
27        (diff0, diff1) = (self.__derivative[0], self.__derivative[1])
28        alpha0 = lambda x: (1 + 2*((x-x0) / (x1-x0))) * ((x-x1) / (x0-x1))**2
29        alpha1 = lambda x: (1 + 2*((x-x1) / (x0-x1))) * ((x-x0) / (x1-x0))**2
30        beta0 = lambda x: (x-x0) * ((x-x1) / (x0-x1))**2
31        beta1 = lambda x: (x-x1) * ((x-x0) / (x1-x0))**2
32        H = y0*alpha0(x) + y1*alpha1(x) + diff0*beta0(x) + diff1*beta1(x)
33        return H
34
35
36if __name__ == '__main__':
37    print('\033[H\033[J')
38
39    obj = HERMITE(
40        X = [0, 1],
41        Y = [0, 1],
42        derivative = [-1, -4]
43    )
44
45    f:'function' = lambda x: obj.Interpolate(x)
46
47    nodes:list = [-0.25, 0.25, 0.75, 1.25]
48    predictions:list = list(map(f, nodes))
49
50    print(f'节点: {nodes}\n\n插值: {predictions}')